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Abstract. The concentration profiles of monomers and counterions in star-branched polyelectrolyte micelles 
are calculated through Monte-Carlo simulations, using the simplest freely-jointed chain model. We have 
investigated the onset of different regimes corresponding to the spherical and Manning condensation of 
counterions as a function of the strength of the Coulomb coupling. The Monte-Carlo results are in fair 
agreement with the predictions of Self-Consistent-Field analytical models. We have simulated a real system 
of diblock copolymer micelles of (sodium-polystyrene-sulfonate)(NaPSS)- (polyethylene- propylene) (PEP) 
with f=54 hydrophilic branches of A = 251 monomers at room temperature in salt-free solution and 
compared the calculated form factor with our neutron-scattering data. 

PACS. PACS-key 82.70.-y - PACS-key 61.20.-p - PACS-key 82.35.Rs 



1 Introduction 

Charged diblock copolymers with a long and fully charged 
hydrophilic part and a short hydrophobic head are attract- 
ing a great attention. In water, their hydrophobic cores 
tend to aggregate and they form colloidal solutions of 
star-branched spherical micelles. These physical systems 
have many applications in surface adhesion, stabilization 
of colloidal solutions and are also used in drug delivery 
processes. Due to the presence of long-range Coulomb in- 
teractions, simple scaling arguments are generally unsuffi- 
cient to describe their behavior. Quantitative predictions 
concerning the concentration profiles of monomers and 
counterions have been obtained through Self-Consistent- 
Field (SCF) models Numerical simulations can 
also provide useful information. The Molecular Dynam- 
ics (MD) approach has been widely used in the study of 
\ neutral star-branched micelles and charged linear poly- 
electrolytes f||6). Recently Jusufi, Likos and Lowen (JLL) 
applied it to star-branched polyelectrolyte micelles. 
They studied, at relatively low linear charge densities of 
the chains, the monomer and counterion profiles of one 
micelle and the effective interaction between two micelles. 
They also proposed analytical approximations to be com- 
pared with the results of MD simulations. 

We focus here our attention on the problem of one 
isolated micelle, in a different strongly-charged regime, 
which is relevant for some dilute solutions of (sodium- 
polystyrene-sulfonate)(NaPSS)- (polyethylene- propylene) 
(PEP) that have been investigated through neutron scat- 
tering in our laboratory B. 



We use the Monte-Carlo method which, through a stochas- 
tic approach, should provide the same information as the 
Molecular Dynamics. Our final goal is to reach real-system 
sizes like stars with / = 54 arms of N = 251 monomers 
(27108 particles), corresponding to experimental condi- 
tions and compare with our neutron-scattering data, with- 
out any adjustable parameter. As a prelude, we first present 
a systematic investigation of smaller micelles within a 
more general theoretical framework and compare our nu- 
merical Monte-Carlo results with the predictions of scaling 
theories and analytical approximations. 

After the description of our model and Monte-Carlo 
method (section 2), we study in section 3 the evolution of 
the structure of a micelle with / = 27 arms of N = 130 
monomers as a function of the following dimensionless pa- 
rameter measuring the strength of the Coulomb interac- 
tion: 

C = N£ = \zi\\z m \l B /a (1) 



where 



l B = (e 2 /4iree )(l/k B T) 



(2) 
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is the Bjerrum length, £ = |z m |/s/d represents the "Man- 
ning parameter" |L(| which characterizes the linear charge 
density of polymers; a is the distance between consecutive 
charged monomers, T is the temperature, e the dielectric 
constant of the solvent and z m , zi represent the valence of 
monomers and counterions respectively. Throughout our 
study z = —Zm — Zi > and all monomers are charged. 

The Bjerrum length is Ib ~ 0.71 nm in water at room 
temperature and it can be varied through the temperature 
and dielectric constant. The monomer length is usually 
of the order of a fraction of a nanometer. £ is approx- 
imately 3 for sodium polystyrene sulfonate in water. It 
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can be increased by choosing other solvents with lower di- 
electric constants and/or by taking polyvalent monomers 
and counterions (z > 1). Higher values of £ are thus eas- 
ily accessible. Although, in our fully charged model small 
values of £, much less than 1, are at present experimen- 
tally unrealistic, it is nevertheless interesting to extend 
our theoretical investigations down to that range. 

When the thermal energy is much larger than the Cou- 
lomb interaction (£ << 1), most counterions are outside 
the micelle corona. This corresponds to the "unscreened 
regime" Q. Above some threshold value ( sp h ~ 1//, a 
spherical condensation of counterions inside and in the 
vicinity of the corona is expected p],^,p|,[T]]. In this so- 
called "osmotic regime" the spherically condensed 
counterions lower the interaction between branches. Both 
unscreened and osmotic regime have been extensively stud- 
ied through mean-field analytical models [|I}||,||,|l2]]. Fi- 
nally, when the "Manning criterion" (To) £m ~ l/^i, i.e 
Cm 1, is reached, a "cylindrical condensation" of part 
of the counterions around the branches is superimposed to 
the spherical condensation. Scaling and SCF theories jl], 
||,||,[l3) have predicted a linear extension of the branches 
with corona radius R c proportional to N and monomer 
density varying roughly as 1/r 2 , where r represents the 
radial distance to the center of the micelle. This has jus- 
tified a "rod-like" picture of the branches. This picture 
has been confirmed by clear experimental evidence of the 
stretching of the chains jL4|. For highly charged polymers, 
a Poisson-Boltzmann analysis jl6| has been proposed to 
interpret the distribution of counterions around the rod- 
like branches in terms of a ionic condensation |Io| , p^| . 

In the three latter regimes, we have compared the re- 
sults of our simulations to the expectations from approxi- 
mate analytical theories. The simulation of a real system 
of (sodium-polystyrene-sulfonate) (NaPSS)- (polyethylene- 
propylene)(PEP) with / = 54 branches of N = 251 mono- 
mers, is described in section 4, with a direct comparison 
to our neutron-scattering data M . 



2 The Model and Monte-Carlo method 
2.1 The model 

We use the simple "freely-jointed chain" (or "pearl neck- 
lace") model. Each monomer is represented by a hard 
core of diameter a m = a. The angle between two suc- 
cessive monomer-monomer segments is free. We consider 
that the hydrophilic chain is fully ionized in water and 
each monomer carries a charge (z m e). This is relevant for 
sodium- polystyrene- sulfonate, in our experimental con- 
ditions, where the degree of ionization is generally of order 
0.95. The counterions are modeled by hard spheres of di- 
ameter Ui = a/2, they carry a charge (^e). For sodium- 
polystyrene- sulfonate Zi — —z rn = 1. 

One isolated micelle and the counterions are treated 
in the "primitive cell model" i.e. they are confined inside 
a hard-wall cube of edge A, surrounding the center of the 
micelle. The dielectric constant is assumed to be the same 



inside and outside the cubic cell. This approach is rele- 
vant for low concentrations with respect to the overlap 
concentration c* of micelles. The weak dependence of our 
physical results on the size A of this cubic box (provided 
it is larger than the diameter of the micelle) means that 
a more sophisticated treatment of Coulomb interactions 
(e.g. periodic boundary conditions and Ewald sums) is 
not required here. 

The hydrophobic core is modeled by a hard sphere 
of radius R CO re of order 5 to 20a. The hydrophilic parts 
are attached at fixed points chosen on this sphere. Two 
distributions of these fixed points have been considered: 

— a random distribution 

— a regularized distribution with approximate equidis- 
tance between first- neighbor pairs. 

The results obtained with these different distributions are 
practically identical and it is thus not useful to perform 
mean values of Monte-Carlo results over different realiza- 
tions of random distributions of points on the sphere. 

All distances are expressed in units of the monomer 
length a and the strength of the Coulomb interaction com- 
pared to ksT is measured by the dimensionless parameter 
£ defined in equation (1). 



2.2 Monte Carlo moves 

For a counterion, an elementary move is chosen at random 
in a cube of edge Sx rs 10a centered around it. 

For polymers we use the "pivot algorithm" , which has 
been shown to satisfy ergodicity |l7| . A monomer P is cho- 
sen at random, and a unit vector u is taken randomly in 
some solid angle 5f2 centered about the axis defined by P 
and the next monomer (P+ 1). A rigid rotation of random 
angle around u is applied to the part of the polymer going 
from (P + 1) to its end. However, beyond Manning con- 
densation (£ > 1), such a pivot algorithm involving only 
the monomers becomes inefficient since many counterions 
are localized around each polymer and the cost in energy 
to move a group of monomers away from these counterions 
prohibits any Monte-Carlo move. 

To circumvent this problem, we have used the "par- 
tially clothed" pivot algorithm proposed by Gordon and 
Valleau JlJ] . Having chosen a pivot P and the correspond- 
ing end part of the polymer, we look for all counterions at 
a distance less than d of this polymer part (the distance 
from the polymer is defined as the distance to the closest 
monomer). Each of these counterions is rotated in block 
with the part of the polymer, with a probability p < 1. The 
detailed-balance equation is satisfied if, in the Metropolis 
scheme, we take as the probability of acceptance of the 
Monte-Carlo move [[jl] 

P acc = mm{l, (i-p)^««-"-w cxp [-(U new - U old )/k B T}}, 

(3) 

where n new and n id are the number of counterions at a 
distance less than d after and before the temptative rota- 
tion. U new and U id are the Coulomb energies after and 
before the trial move. 



M. Roger et al.: Monte-Carlo simulations of star-branched polyelectrolytes 



3 




-1.0 1 1 1 1 1 1 1 1 1 

500 1000 1500 2000 

Monte-Carlo steps per particle 

Fig. 1. Coulomb energy as a function of the number of steps 
at £ = 2 . Full black line: starting from fully stretched poly- 
mers. Dashed line: starting from the configuration of neutral 
polymers. Full gray line: starting from the equilibrium config- 
uration at C = 1. Here N = 130, / = 27, A = 320a. Beyond 
1000 steps the differences between the three curves are of the 
order of the fluctuations. 
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Fig. 2. Number of accepted "cluster-moves" per monomer as a 
function of the monomer label (from 1 to N, starting from the 
core). Full line: "clothed pivot rotations"; dashed line: "clothed 
translations" (see text). 



In practice, we have observed that "clothed pivots" 
have a reasonable probability to be accepted in the Metro- 
polis scheme if p > 0.9. Since we did not observe very 
important differences in the convergence to equilibrium 
for 0.9 < p < 1, we chose the value p = 1 which makes the 
algorithm simpler and the calculation shorter. (The choice 
p = 1 has also been made by Lobaskin and Linse jyi 
in similar "cluster moves" of macro-ions surrounded by 
a shell of small counterions). The maximum ion-polymer 
distance for clothed pivots has been optimized to d = 4a. 

To facilitate stretching and contraction of the polymer, 
we have also applied some small "clothed block transla- 
tions" of parts of polymer in the following way. A monomer 
P is chosen at random. The next monomer (P + 1) is 
moved from rvp+ji to r' P , 1 within a solid angle 8fl cen- 
tered around (r^ P+1 - j — rp). The rest of the polymer from 
P to its end is translated in block by a vector (r^ p+1 * j — 

r (p+i))- Such moves are also combined with "cluster moves' 
of counterions at distance less than d. Compared to piv- 
ots they have a higher probability to be accepted and a 
suitable combination of these both types of cluster-moves 
accelerates the convergence. 

Our main criterion for equilibrium is the convergence 
of the Coulomb energy to the same stationary value start- 
ing from two drastically different initial configurations: 

— Fully extended polymers with counterions distributed 
uniformly in the volume of the cell. 

— Contracted star with polymers in the equilibrium con- 
figuration of uncharged polymer and counterions dis- 
tributed uniformly in the volume of the cell. 



As a typical example, Fig. 1 illustrates the convergence 
of the Coulomb energy, starting from these two extreme 
configurations for a micelle with 27 branches of 130 mono- 
mers (2 x 3510 particles) at £ = 2. The equilibrium was 
reached after 10 3 steps per particle, representing 25 hours 
of monoprocessor time on a Compaq RS232 parallel com- 
puter. The relevant statistical properties (distribution of 
monomers and counterions, form factors etc..) were mea- 
sured over 10 3 further steps. A "cluster move" of a poly- 
mer part was attempted every 3 iterations (the two oth- 
ers being devoted to single counterion moves). The cluster 
moves were either pivots (3/4 of them) or translations (1/4 
of them) as previously described. During the measure of 
statistical averages, each counterion was moved approxi- 
mately 1300 times, about 80 pivots and 200 translations 
per monomer were accepted. 

The number of accepted pivots as a function of the 
position of the pivoting point P is represented in Fig. 2. 
Close to the core of the micelle, the number of accepted 
moves is low (of the order of 10 per monomer). Neverthe- 
less, it grows roughly exponentially as a function of the 
position of monomer P and statistics becomes acceptable 
beyond a few tens of monomers. 

The size Sx of the cubic box for elementary counterion 
displacements was 10a, and the elementary solid angle 8f2 
for pivot rotations was fixed to 27rl0 -2 steradians. We 
have also tried some simulation with 5f2 increasing with 
the distance to the core of the micelle, without any sub- 
stantial improvement of the rate of convergence to equi- 
librium. 
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3 General results 

3.1 The theoretical model studied 

We present here some theoretical results concerning the 
variations of the monomer and counterion profiles as a 
function of the dimensionless parameter £ [cf. equation 
(1)] measuring the strength of the Coulomb interaction 
with respect to fceX '. 

We consider two micelles with numbers of monomers 
and branches smaller than these corresponding to the ex- 
perimental system studied in the next section: 

(i) N = 130 and / = 6 

(ii) N = 130 and / = 27 

The radius of the core R CO re is chosen to be small com- 
pared to the radius of the micelle corona and has thus a 
little influence on the corona. We take in all cases R CO re = 
5a. The influence of the boundary conditions was investi- 
gated by comparing the results obtained with two different 
sizes of the cubic cell: A = 320a and A = 640a. 

We vary the dimensionless parameter £ keeping the ra- 
tios A = A/ a, Rcore — Rcore/a and &i — Oij a constant. 
We start from a low value Co = 0.005 and increase it by 
steps. At each new step Cn+i we take as initial state of the 
Monte-Carlo calculation the equilibrium state obtained at 
the preceding one £ n . The system is thus progressively 
"annealed" , which provides a faster convergence to equi- 
librium at each step. As a typical test for equilibrium we 
compare in Fig. 1 the energy obtained in that way at £ = 2 
starting from the equilibrium configuration at C = 1 to 
the results obtained with far-from-equilibrium arbitrary 
initial configurations. After 1000 iterations the differences 
between all curves are of the order of the fluctuations. 

We calculate the mean values of the radial distribu- 
tions of monomers and counterions. The deviations of these 
distributions from the monomer distribution correspond- 
ing to fully extended arms are better emphasized by con- 
sidering, in terms of f = r/a, the dimensionless quantities: 

P{m,i}{r) = — r 2 a 3 c { „ M} (f) (4) 

where c/ m n(f) represent the local volume densities. 
For a micelle with fully extended arms, p m should be 1 
for R core < r < R core + N and elsewhere. 



3.2 Monomer profile as a function of ( 

For N — 130, f — 27 and A — 320, typical monomer and 
counterion profiles at represented in Figs. 3a, 3b and 3c at 
weak coupling (£ = 0.02), intermediate coupling (( = 0.2) 
and strong coupling (£ = 3) respectively. Corresponding 
snapshots of the micelle are shown in Fig. 4. 

The variations of the radius of gyration R g and of the 
corona radius R c as a function of £ are shown in Fig. 5 
(circles). The corona radius R c is defined here as the ab- 
scissa of the maximum of the free-end distribution, which 



is approximately gaussian (see Fig. 6). It also roughly co- 
incides with the rightmost inflexion point in the concen- 
tration profile p m (r). 

We have studied the influence of boundary conditions 
by doubling the size A of the cubic cell (crosses). R c and 
R g are very weakly dependent on boundary conditions. 
Relative difference of the order of a few percent at most 
are observed at maximum stretching (£ of order 1). We 
conclude that it is not useful to investigate more sophis- 
ticated treatments of Coulomb interaction like these cor- 
responding to periodic boundary conditions with Ewald 
sums. 

The influence of the number of arms is seen in Fig. 5 
by comparing our results with a small number of arms 
/ = 6 (diamonds) to the results with / = 27 (circles). It 
is only important at weak coupling. 

The Figs. 3, 4 and 5 show clearly, as a function of £, 
three different regimes. These three regimes are governed 
by changes in the behavior of the counterions, which are 
expected from mean- field theories |J[§|,|h|[11]]- 

3.3 The "unscreened regime". 

At weak coupling £ << 1, the contribution of the coun- 
terions to the entropy dominates the electrostatic energy. 
The local volume density a of counterions is roughly con- 
stant in the cell: Ci w fN/A 3 and from equation (4) 
Pi w AnNf 2 1 A 3 . If A is larger than R c , the fraction of 
counterions in the micelle corona is small and the electro- 
static interactions between branches are weakly screened. 

In the unscreened limit, the self-consistent field argu- 
ments proposed a long time ago for linear polyelectrolytes 
by de Gennes et al. f P- 3| have been straightforwardly gener- 
alized to star-branched polyelectrolytes 0,^]. Let us con- 
sider a chain whose first unit is placed at the origin. The ra- 
th unit is at an average distance r(n) and feels two forces: 

— In a simple gaussian model, there is an elastic force 
which tends to contract the chain: 



F,, 



3k B T d 2 r 
a 2 dn 2 



(5) 



— If we assume equal stretching of the / chains of a star, 
i.e. all chain ends are at the same distance R c from the 
origin, the spherical region at distance smaller than 
r{n) contains a mean charge Q = nfz m e. Applying 
Gauss theorem, there is an electrostatic force: 



F n = 



f{z m e) 



1 



47reeo r(n) 2 



(6) 



Writing the equilibrium between the two forces, one ob- 
tains the following dimensionless equation: 



d 2 u 1 x 
dx 2 3 u 2 







(7) 



with x = n/N and u — r/(Nd) where the length scale d 
is: 



d = 



(8) 
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Fig. 3. Density profiles of monomers (black) and counterions 

(grey) in the three different regimes: a) uncondensed regime Fi S- 4 - Snapshots of the micelle in equilibrium corresponding 
(C = 0.02); b) spherically condensed regime (C = 0.2); c) to the three regimes considered in Fig. 3 respectively. The three 
Manning-condensed regime (C = 3). Here N = 130, f = 27 snaphots are represented at the same scale, 
and A = 320. 
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Fig. 5. Variation of the corona radius R c and radius of gyra- 
tion R g as a function of the dimensionless parameter f for: (i) 
/ = 27, A = 320 (circles); (ii) / = 27, A = 640 (crosses); (iii) 
/ = 6, A — 320 (diamonds). N = 130 in all cases. The Manning 
threshold is indicated by the dash-dotted vertical line £ « 1. 
The dashed and dotted vertical lines indicate different thresh- 
old values for spherical condensation corresponding to (i) , and 
(iii) respectively . The horizontal arrow on the left indicate the 
limiting values of the corona radius corresponding to neutral 
stars with / = 27 arms of N = 130 monomers. 



Fig. 6. Distribution of branch ends for £ = 0.02. The Monte- 
Carlo results (full line) are fit by a gaussian (dashed line): 
l/[v / 2T(?)] exp[-(f- J R c ) 2 /2? 2 ] with R c = 66.33 and a = 5.31. 




Fig. 7. Monomer profiles in the "unscreened" regime. Full 
lines: N = 130 and / = 27 for three (-values, < = 0.005 (lowest 
maximum), £ = 0.01 and £ = 0.02 (highest maximum). Dotted 
lines: N = 130 and / = 6 for < = 0.01 and C = 0.02. The gray 
dashed line indicates the solution of Eq. 7. 



M. Roger et al.: Monte-Carlo simulations of star-branched polyelectrolytes 



7 



The boundary conditions are u(0) = and du/dx = at 
X = 1. Taking / = 1 in the previous algebra is nothing 
else than De Gennes et al. equations |l3} |. 

The center to end distance of the chains is thus pro- 
portional to Nd: 



1.0 



Rc ~ f lf \ 1/3 Na 



(9) 



(we take Zi = —z m ). 

These simple scaling arguments lead to a stretching 
of the chains which is proportional to the number TV of 
monomers. They have justified a so called "rod-like" pic- 
ture of the branches. 

In Fig. 5, the slope d log R c /d log ( is smaller than 
1/3, the value expected through the previous scaling ar- 
guments. The dependence of R c on / is also weaker than 
predicted by equation (9). At very weak coupling, we ob- 
tain (see Fig. 5) R c (f — 27)/ R c (f = 6) « 1.35 compared 
to the expected scaling: (27/6) x / 3 = 1.65. The explanation 
is that when £ is low enough to ensure that very few coun- 
terions are inside the corona, the excluded volume effects 
(which are not taken into account in the previous scal- 
ing arguments) become significant. The value of R c has a 
lower limit which is given by the scaling law for neutral 
stars @: 

R neutral _ ^1/5^3/5 (10) 

This limiting value for N — 130 and / = 27, within our 
pearl-necklace model is indicated in Fig. 5. Our lowest 
value for R c shows an inflexion towards this asymptote. 

Close to the core, our monomer profiles can be roughly 
scaled on a universal curve by plotting p m (R c — R core )/Na 
as a function of [r—R core ) / (R c — R core ) (i.e. the integral of 
each curve is normalized to 1) and compared to the solu- 
tion of Eq. 7 given in Fig. 2 of Ref. (see Fig. 7). Close to 
the core, our Monte-Carlo results agree qualitatively with 
the solution of equation (7) which gives a roughly constant 
profile of p m (i.e. a 3 c m (r) ss f~ 2 ). At R c , equation (7) 
leads to some unphysical divergence, which is smoothed 
in more refined SCF models taking into account the fluc- 
tuations of the chain ends jl],|j| . Nevertheless a substantial 
increase of p m near R c remains in these SCF models, as 
well as in our simulation. 



3.4 The "osmotic regime" with spherical condensation 
of counterions 

Increasing C we expect a "spherical condensation" of a 
part of the counterions inside and in the vicinity of the 
corona. This occurs when the Coulomb energy of a coun- 
terion in the field created by the charged monomers: 

Nf) Zl e 2 /{4iree Rc) 

at the limit of the corona radius is of the order of kT, i.e. 
for 

Csph 



C 

5> 
o 

3= 
CD 
O 

o 
o 



CO 

o 




0.0 



Fig. 8. Osmotic coefficient 0. The open diamonds and circles 
represent the Monte-Carlo results for / = 6 and / = 27 re- 
spectively. The full and dashed line correspond to the Poisson- 
Boltzmann approximation. In both cases A = 320 and N = 
130. 



(at this level of approximation we take R c k Na). The 
condensation of ions near the surface of an impenetra- 
ble charged sphere 1 1 1 has been formulated in an analo- 
gous way to the cylindrical condensation of ions around 
a charged rod studied by Manning Jic[| . The more com- 
plex problem of counterion condensation in a penetrable 
micelle has been studied through Self-Consistent Field 
methods @,|,|. It can also be formulated in a simpler 
way through a Poisson-Boltzmann analysis. The Poisson- 
Boltzmann equation can be solved numerically for a pen- 
etrable sphere of radius R c with charge density 



/ 



c m (r) = fN/(4irR c r 2 ) 

inside a spherical cell of volume A 3 . The Fig. 8 compares 
the Monte-Carlo to the Poisson-Boltzmann results for the 
osmotic coefficient c/>. The osmotic coefficient is defined as 
the ratio of the osmotic pressure to its value corresponding 
to a uniform repartition of the ions in the cell. The osmotic 
pressure is equal to ksT times the local ionic density at 
the edge of the cell. 

There is a clear change of slope of c/>(C) at a point which 
can be defined as the crossover to the condensed regime. 
For different numbers of arms and cell sizes, we deduce 
the crossover to the spherically condensed regime at the 
following values Csph'- 

- Csph « 0.15 for / = 6, A = 320 

- Csph ~ 0.03 for / = 27, A = 320 

- Csph w 0.08 for / = 27, A = 640 

(N = 130 in the three cases). At fixed A, Csph varies in 
!//• 
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Fig. 10. Variation of the corona radius near and in the Man- 
ning condensed regime. The stars represents the results of JLL 
[7] for N = 50 and / = 18 at charge fraction a = 1/6, 1/4, 1/3, 
the triangle corresponds to their results for N = 150, / = 10 
and a — 1/3, the abscissa is f = a£. In JLL [7], all simula- 
tion are at constant ( — 3. The other symbols are our results 
reported in Fig. 5 with the same conventions (in our study a 
is fixed to 1 and £ = f ' is varied) . The dashed arrows indicates 
a mapping of our model at £ = 3 and a = 1, in the strongly 
condensed Manning regime, onto an effective condensed sys- 
tem with charge fraction a e ff = 1/3. It agrees nicely with the 
real system at a = 1/3 studied by JLL at the same f = 3. 



In Fig. 9, we show that at the thresholds defined above, 
there is a drastic change of the counterion profile. Below 
Csph, the local density Ci of counterions is roughly constant 
Ci « fN/A 3 leading to a quadratic behavior of pi. Above 
the threshold, the curves pi show an inflexion point at 
r ps R c ja. Since in that £ range R c is roughly constant, 
the abscissa of this inflexion point remains constant. 

For C >> 1//, most of the counterions are inside the 
micelle corona (see Fig. 9). This regime has often been 
called "osmotic" by reference to the osmotic pressure of 
counterions inside the corona |^,||. 

The range 1// << ( < 1 corresponds to the largest 
stretching of the star with a flat maximum in the curve 
R C {Q (see Fig. 5). 



Fig. 9. Radial distribution of counterions (N=130) as a func- 
tion of ( for different numbers of arms: (a) / = 6 and A = 320, 
(b) / = 27 and A = 320, and different cell sizes: (c) / = 27 
and A — 640. The arrows indicate a common inflexion point at 
t ~ He in the pure spherically-condensed regime. The dashed 
parabola indicate the equi-partition of all counterions in the 
volume of the cell in the uncondensed regime. 



3.5 The "Manning condensed" regime 

For larger £ values, beyond the "radial condensation" of 
counterions in the micelle corona, we expect the superposi- 
tion of a "cylindrical condensation" of counterions around 
each branch of the micelle. This phenomenon is clearly 
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Fig. 11. Monomer density profiles for ^ — 2, N = 130 and 
/ = 6 (dotted line), / = 27 (dashed line), and / = 54 (full 
black line). They are compared to the density profile of one 
branch attached to a fixed center point with R c — (grey 
line). Since R c — 5a for all other curves this last grey line has 
been translated by 5 units in x direction. The dot-dashed line 
indicate the profile corresponding to fully extended branches. 
Only a slight increase of the corona radius is observed when 
the number of branches is varied by a factor 10 (from / = 6 to 
/ = 54) 

seen in Fig. 3c and Fig. 4c. The Manning criterion for this 
"cylindrical condensation" is |io| , pd| : 

6/ = 1/1*1 

which corresponds to the threshold value (m = 1< 

The Fig. 10 represents an enlargement of the last decade 
in Fig. 5. The corona radius R c presents a flat maximum 
at 0.5 < ( < 1 and decreases at ( > 1. At this stage, it 
is interesting to compare our approach with the Molecu- 
lar Dynamics study by JLL |t],|| where £ has been fixed 
to 3, but where only one monomer out of three, four or 
six is charged. The "charge fraction" a is thus varied 
from a = 1/3 to a = 1/6 whereas it is a = 1 in our 
model. The relevant Manning parameter which character- 
izes the linear charge density of the polymer is £ = z m Zs/6, 
where b = a/a is the distance between consecutive charged 
monomers pjj. In Fig. 10 we have plotted some of JLL re- 
sults, as a function of £' = a(. The fact that from a = 1/6 
to a = 1/3 their corona radius increases indicates that 
they have not reached the Manning condensed regime at 
a = 1/6 and a — 1/4. They are just reaching it near 
a = 1/3. 

At £ >> 1//, (near the broad maximum of R c and in 
the Manning condensed regime) , the monomer profiles are 
only very weakly dependent on /. In Fig. 11, we compare 
the monomer profiles for different values of the number 
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Fig. 12. Variation of the gyration radius R g as a function of 
the number of branches / for £ = 2 and N = 130. Only an 
increase of a few percent of R g is observed when the number 
of branches is varied by a factor 10 



of branches / = 6, 27, and 54. The number of monomer 
per branch is N = 130, the radius of the hydrophobic 
core is R core = 5 and the size of the cell is A = 320. The 
Fig. 12 represents the variation of the gyration radius as 
a function of the number of branches /. It increases only 
by a few percent when the number of branches vary by a 
factor of 10. The MD results of JLL (see Tables I and II 
in Ref. ||), also show a weak dependence of R c on /. 

The Fig. 13 represents a blow up over the last 80 mono- 
mers of one branch in Fig. 4c. At a length scale of the 
order of 10 monomers, a local "rod-like" picture is rele- 
vant. To illustrate the phase space explored by a given 
branch, we have registered 100 snapshots of the same 
branch (Fig. 14), at equal time intervals in thermody- 
namic equilibrium. The phase space explored corresponds 
roughly to a cone of solid angle 4.7r// steradians which 
shows that each branch occupies a fraction 1/f of the 
sphere and rarely "mixes" with neighboring branches. An 
"urchin-like" picture |2^| is thus relevant. 

The last decade in Fig. 5 and its enlargement in Fig. 10 
show that, in our fully charged model (a = 1), beyond 
the Manning threshold, the corona radius R c decreases in 
£-1/3 rp n j g - g q U jj. e a new resu it. Both this exponent and 

the independence of R c on / can be interpreted with the 
following simple arguments. 

For £ > 1, there is a cylindrical condensation of a frac- 
tion of counterfoils along the branches. The model pro- 
posed by Manning [jl(|, considers that a fraction (1 — l/£) 
of the counterfoils condenses around the branches in such 
a way that the monomers get an effective charge a e ffz m 
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Fig. 13. A blow-up over the last 80 monomers of one branch 
in Fig. 3c. Large spheres are fully charged monomer. Small 
spheres with diameter 1/2 of the large ones are counterions 



with an effective charge fraction 

a eff = 1/C- 

The remaining fraction 1/C of counterions are only spher- 
ically condensed and distributed between the branches. 
We assume that / is sufficiently high (C >> 1//) to ne- 
glect the fraction of "free" counterions outside the corona. 
Taking a spherical distribution in 1 jr 2 , we obtain the fol- 
lowing effective ion density: 



1 



(r) = 7 



fN 



C A-Ka 2 R c f 2 ' 

The corresponding "effective" screening length is 



if) 



/t 



Na 



Rc 



-1/2 



Since we restrict our study to \zi\ = \z rn \, it is inde- 
pendent of £. The "effective" screening length roughly 
keeps the same value it reaches at the Manning conden- 
sation threshold Q — 1. At a radial position r, the mean 
transverse distance between two neighboring branches is 
roughly 4af/ -1 / 2 . It exceeds the screening length by a 
factor larger than 4. Hence we can neglect the electro- 
static interactions between branches. Since the polymers 
are relatively extended, the excluded volume interactions 
between branches are also weak. This is the reason why 
the monomer profiles are so weakly dependent on /. 

To justify the (—1/3) exponent, we consider one branch 
in the familiar electrostatic "blob" picture jl3). Let us rep- 
resent a branch as a succession of N/m blobs of size I, 




Fig. 14. Snapshots of the same branch at equal time intervals 
during measurement of thermodynamic quantities. The portion 
of space explored corresponds to a cone of solid angle 4n/f 
steradians. The sphere at the top represents the hydrophobic 



each containing m monomers. As far as the size of a blob 
is smaller than the screening length k ~^ , the electrostatic 
energy of a blob is 

E c = (a e ffz m em) 2 /(4iree l). 

Its elastic energy from entropic origin is 

E e = (3/2)k B Tl 2 /(ma 2 ). 

Minimizing the sum E c + E e with respect to m, we obtain: 



(I /ma) 



4 (a e ^) 2 z 2 m l B 
3 a 



1/3 



4 | Z-m | 

3N" 



-1/3 



and the length of the branch is 

R c = (N/m)l = Na 



4 \z m \ 
3 \zi\ 



-1/3 



\zi\, R c is proportional 



Since in the present study \z m 
toC 1/3 {§. 

At high Schiessel and Pincus predicted a col- 
lapse of the polymers with the following arguments. When 
the Coulomb energy dominates over the entropic contribu- 
tion, each counterion is practically fixed to a monomer and 
forms an electric dipole. The attraction of dipolcs should 
lead to a transition towards a collapse of the star. Up to 
C = 10 our results follow a C -1 / 3 power law, and we do 
not observe any precursor of this effect. 

The Manning picture allows us to make some approxi- 
mate mapping of our model at a = 1 , ( = 3 onto the model 
of JLL at a = 1/3 and the same C = 3. In our model, at 
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Fig. 15. Density of counterions at a radial distance f from 
the micelle center and a transverse distance from the nearest 
branch d. The results of the Monte-Carlo simulation (symbols) 
are compared to those of a Poisson-Boltzmann approximation 
(continuous lines), (a) Circles and dotted line: f = 20; (b) stars 
and dashed line: f = 40; (c) diamonds and full line: f = 70. 
Here, C = 3, N = 130, f = 27 and A = 320. 

( = 3, due to Manning condensation, (1 — 1/C) = 2/3 of the 
counterions condense on the monomers. The remaining ef- 
fective charge fraction of monomers is thus a e // = 1/3 and 
the other one third of the counterions is only spherically 
condensed inside the corona. This effective model is thus 
close to the real situation studied by JLL at a = 1/3 and 
C = 3. The arrow in Fig. 10 shows a very nice agreement 
of the corona radii corresponding to this mapping. This is 
a strong evidence of the relevance of Manning's picture. 

Up to now, radial distributions of counterions were 
only considered. In the Manning regime, it is of prime in- 
terest to study the cylindrical distribution of counterions 
around one branch. Defining the distance d of a counte- 
rion to the closest branch as the distance to the closest 
monomer we have measured the mean value of the ra- 
dial density n(r,d), n(r,d)SrSd representing the number 
of counterions at distance between d and (d + Sd) from a 
branch and at distance between r and (r + Sr) from the 
center of the micelle. 

The Fig. 15 represents the density n(f, d) of counteri- 
ons at distance d = da from the closest branch and at a 
distance fa from the center of the micelle for three values 
of f: close to the core (f = 20), in the middle (f = 40) and 
close to the edge of the corona (f = 70). We chose £ = 3, 
well above the Manning threshold, N = 130, / = 27 and 
A = 320. 

It is interesting to compare these results to those ex- 
pected from a "rod-like" picture of the arms and a Poisson- 
Boltzmann approach for the distribution of counterions. 



Fig. 16. Radial density of counterions. The Monte-Carlo re- 
sults are compared to the Poisson-Boltzmann approximation 
(dashed line). Here, C = 3, N = 130, / = 27 and A = 320. 




d 

Fig. 17. Number Ui(d) of counterion per monomer at distance 
less than d from the closest branch, for different values of the 
dimensionless parameter £. Here N = 130, / = 27 and A — 
320. The change of curvature between £ = 0.5 and £ = 1.5 
corresponds to Manning condensation. The Manning threshold 
£ = 1/zi is reached at ( f» 1. A logarithmic dependence of 
v on d is expected at £ = 1/zi, from a Poisson-Boltzmann 
approximation in a "rod-like" picture. For each value of £ in the 
Manning condensed regime, the arrow indicates the distance at 
which the "effective charge parameter" a e f / = 1 — i>i leads to 

Zeff = 1 / z i- 



12 



M. Roger et al.: Monte-Carlo simulations of star-branched polyelectrolytes 
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Fig. 18. Dependence of Vi(d) on the number of arms / and 
boundary conditions (size A of the cubic cell). Full lines: TV = 
130, f = 27 and A = 320; dashed lines: N = 130, / = 27 and 
A = 640 ; dash-dotted lines: N = 130, / = 6 and A = 640. 

Using the "FIDISOL/CADSOL" package developed at 
Rechenzentrum, Universitat Karlsruhe for the numer- 
ical solution of non-linear partial differential equations, we 
have solved the Poisson-Boltzmann equation for a rod in a 
cone of solid angle (An/f). The radius of the rod is chosen 
to 0.75a, which corresponds to the sum of the hard-core 
radius of monomer and counter ion. The height of the cone 
is equal to the radius of a sphere of volume A 3 . The distri- 
bution of charges in the rod corresponds to the monomer 
profile obtained in Fig. 3c through Monte-Carlo simula- 
tion. (In practice, we fit this distribution to an approxi- 
mation by a rational function). Using a two-dimensional 
80 x 120 grid, we obtained an accuracy better than 10~ 4 on 
almost every point. The calculation took only one minute 
of monoprocessor time. There is a nice agreement except 
in the vicinity of the corona radius. 

Fig. 16 compares the overall density profile Pi(f) ob- 
tained through the Poisson-Boltzmann approach to the 
Monte-Carlo results of Fig. 3c. 

An interesting quantitative way to characterize the 
Manning threshold is to plot in log-linear coordinates the 
number Vi{d) of counterions at distance less than d from 
the closest monomer (i.e the integral over d' of n(r,d') 
between 3/2 and d) fllqj . Within the Poisson-Boltzmann 
approximation, z^(d) is expected to follow a logarithmic 
law at the Manning threshold |Tq] . As shown in Fig. 17, 
there is a change of curvature of the lines between £ = 0.5 
and ( — 1.5. A straight line (logarithmic behavior) is ob- 
tained at C = 1- 

The Fig. 18 shows that, in the Manning condensed 
regime, Vi(d) is practically independent of the boundary 
conditions: the results obtained after doubling the size of 
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Fig. 19. Density profiles of monomers (black line) and counte- 
rions (grey line) for a micelle with 54 branches of 251 monomers 
at C = 2.8 with R core = 18a and A = 560. 

the cubic cell (A = 640) are practically identical to those 
obtained with A = 320. Vi{d) depends on the number of 
arms but the threshold value Qm ~ 1 corresponding to 
Manning condensation is unchanged. 



4 Simulation of an experimental system 

4.1 Description of the experimental system 

We used asymmetric diblock copolymers made of fully 
dcuterated sodium poly- (styrene-sulfonate) (NaPSSd) and 
of (polyethylene-propylene) (PEP) |9|. These materials 
were synthetized by anionic polymerization. The degrees 
of polymerization were measured and found to be 251 for 
(NaPSSd) moiety and 52 for PEP moiety. The aggregation 
number, i.e the number of branches per micelle, was found 
to be 54 and the core radius is R cor e ~ 4.5 nm w 18a. The 
monomer size for (NaPSSd) moieties is a = 0.25nm. Tak- 
ing Ib = 0.71nm as the Bjerrum length of pure water at 
room temperature, the dimensionless parameter £ is evalu- 
ated to ( = Ib /a = 2.8. The concentration is 0.25c* where 
c* represents the overlap concentration. Our cubic cell was 
chosen such that its volume represents the mean volume 
occupied by one micelle in the solution i.e. A = 560. 

4.2 Monomer and counterion profiles 

The density profiles of monomers and counterions are plot- 
ted in Fig. 19. The corona radius is R c — 174a leading to 
(i? c — R core ) /Na w 0.62 this value has to be compared 
to that obtained in section 3 with N = 130 and / = 54 
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Fig. 20. Monomer-monomer form factor for a micelle with 54 
branches of 251 monomers at £ = 2.8. The experimental results 
(circles) are compared to two Monte-Carlo calculation with dif- 
ferent core size (all other parameters being the same) (i) heavy 
full line: Rcore = 12a; (b) heavy dash-dotted line Rcore = 18a. 
The thin lines indicate analytical asymptotic behaviors. Thin 
full line: the form factors of perfect rods of length R c — R CO re\ 
Thin dashed and dotted lines the form factors for spherical 
distributions of monomers in 1/r 2 between R CO re and R c for 
Rcore = 12a and Rcore = 18a respectively. 



at the same £: (R c — R core )/Na w 0.57. We agree with 
the Molecular Dynamics results of JLL [j7|,|| on smaller 
branch lengths showing that the ratio R c /N is not strictly 
constant but slightly decreases when N increases. The 
monomer and ion profiles are quite similar to those ob- 
served in section 3 for N — 130 and / = 27 (compare 
with Fig. 3c). 



4.3 Form factor, comparison to experimental results 

The monomer-monomer form factor obtained through the 
Monte-Carlo simulation, with the parameters described in 
subsection 4.1 is shown in Fig. 20 (dashed line). It is close 
to the monomer-monomer form factor obtained though 
neutron scattering || except in the range q « ir/R core 
where the size of the core has a large influence on the spec- 
trum. Another Monte-Carlo calculation with a slightly 



smaller core R cor e — 12a, (full line) was performed. It is 
in a better agreement with the neutron scattering results. 

It is interesting to compare the simulated and exper- 
imental results to straightforward analytical approxima- 
tions at low and high q. 

(i) At low q, (q < ir/R core ), the form factor can be rea- 
sonably approximated by that of a 1 /r 2 spherical dis- 



tribution of monomers from r 
zero otherwise: 



Rc 



to r = R c and 



F sph (q) 



Si(qR c ) - Si(qR core ) 

qR c 



(11) 



where Si(z) represents the Sinelntegral function [ [26| 
The thin dashed and dotted lines represent this func- 
tion for Rcore — 12a and R CO re — 18a respectively. 
Up to qa w 0.1 they closely follow the correspond- 
ing Monte-Carlo simulations (heavy full and dashed 
lines). At q « n/R core large oscillations appear in the 
function —Si{qR core ), the first of them being respon- 
sible for a dip observed in the simulation. Close to 
this dip, the form factor is very sensitive to the size of 
the core. We find that our simulation with R cor e = 12 
is closer to the experimental results than that with 
Rcore = 18. The interface between the hydrophobic 
core and the hydrophilic branches might happen to be 
not so smooth and as well defined as assumed in our 
model. For example, some hydrophilic monomers may 
be trapped inside the outer part of the core, leading to 
a lower effective core radius, 
(ii) At high g, the branches can be approximated by rods 
of length L = R c — R cor e with negligible thickness. The 
corresponding form factor is: 



F rod (q) = j I ^Si(qL) 



i(qL/2) 



qL/2 



(12) 



The thin full line represents this function with L = 
R c — Rcore = 156a corresponding to our Monte-Carlo 
estimate. At large q, F rod (q) « [1/ fq(Rc- R cor e)] and 
the position of this line of slope (-1) in logarithmic co- 
ordinates gives a nice estimate of the corona radius. 
In conclusion the corona radius obtained through our 
Monte-Carlo simulation is in agreement with the ex- 
perimental value. 

The peak observed in the experimental data at qa ~ 0.02 
arises from long-range correlations between micelles, mostly 
due to the remaining uncondensed charges. It is a collec- 
tive effect which, of course, is absent in our model. The 
concentration is c = 5 x 10~ 3 g/cm 3 . The corresponding 
q wave vector characterizing the inter-micelle correlations 
is q w 27r[o/V //M] 1/3 where M = 5.7 x 10 4 is the molar 
mass of one diblock copolymer and Af the Avogadro num- 
ber. We obtain qa w 0.016 in reasonable agreement with 
the position of the peak. 



5 Conclusion 

Monte-Carlo simulations appear to be a useful tool for 
understanding the complex structure of highly charged 
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copolymer urchin-like micelles. As the Coulomb coupling 
is increased, three different regimes can be distinguished, 
i) At low coupling, the freely-jointed polyelectrolyte bran- 
ches become stretched due to the monomer-monomer re- 
pulsion and the micelle swells. The monovalent counteri- 
ons are not sensitive to the low miccllar charge and are 
uniformly distributed inside the cell, ii) At intermediate 
coupling, the arms are almost fully stretched and most 
of the countcrions are trapped inside or in close vicinity 
of the micelle. This spherical condensation is due to the 
presence of many, although weakly charged, branches in 
the same micelle, iii) At large coupling, a classical ionic 
condensation appears locally around each highly charged 
branch which adds to the previous one. The effective post- 
condensation monomer charge decreases, the monomer- 
monomer repulsion weakens and the micelle deswells. These 
Monte-Carlo results agree with Poisson-Boltzmann anal- 
yses in different geometries and nicely reproduce small- 
angle neutron scattering data. The influence of intrinsic 
chain rigidity, and the role of ion-ion correlation in pres- 
ence of multivalent counterions will be investigated in fu- 
ture papers. 
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